/**
 * @file FMG_1d.h
 * @brief The FMG for 1 dimension
 * @author XDDDD
 * @version 
 * @date 2021-06-10
 */

#ifndef __PAISLEYPARK__FMG_1D_H__
#define __PAISLEYPARK__FMG_1D_H__

#include "MultigridSolver.h"

class FMG_1d : public MultigridSolver<1> {
public: 
	FMG_1d() : MultigridSolver<1>(){};
	~FMG_1d(){};
	FMG_1d(std::string _RestrictionId, std::string _InerpolationId) : MultigridSolver<1>(_RestrictionId, _InerpolationId){};
	void Solve(VectorFunction<1>& _f, VectorFunction<1>& _g, std::vector<double>& _v,
			   const size_t _N, const size_t _M, const size_t itr_num, 
			   const size_t nu1, const size_t nu2, const double _epsilon) override;
	void VC(const size_t _N, const size_t _M, const size_t nu1, const size_t nu2, const double _epsilon);
};

void FMG_1d::Solve(VectorFunction<1>& _f, VectorFunction<1>& _g, std::vector<double>& _v, 
					 const size_t _N, const size_t _M, const size_t itr_num, const size_t nu1, const size_t nu2, const double _epsilon) {
	bool is_stop = true;
	this->set(_v, _N, _f);
	size_t N = this->__U.get_N();
	size_t pre_N = N;
	for(size_t i = N - 1; i != _M - 1; i--)
		this->set_f(_f, i);
	pre_N = _M;	
	this->VC(pre_N, _M, nu1, nu2, _epsilon);
	if(N > _M)
		is_stop = false;
	while(!is_stop) {
		(*(this->__I_2h_2_h))(this->__U, pre_N);
		pre_N++;
		for(size_t i = 0; i != itr_num; i++)
			this->VC(pre_N, _M, nu1, nu2, _epsilon);
		if(pre_N == N)
			is_stop = true;
	}

	std::vector<double> x(1, 0);
	std::vector<size_t> coord(1, 1);
	size_t pole = 0;
	for(size_t i = 1; i != (1<<_N); i++) {
		x[0] = 0;
		__U(i, _N) += (1.0 - 1.0*i/(1 << _N))*_g(x);
		x[0] = 1;
		__U(i, _N) += 1.0*i/(1 << _N)*_g(x);
	}

};

void FMG_1d::VC(const size_t _N, const size_t _M, const size_t nu1, const size_t nu2, const double _epsilon) {
	bool is_stop = false;
	size_t pre_N = _N;
	size_t l = 1;
	while(!is_stop) {
		for(size_t i = 0; i != nu1; i++) {
				this->smooth(pre_N);
		}
		if(pre_N == _M) {
			is_stop = true;
			break;
		}
		this->set_e(pre_N);
		pre_N -= 1;
		std::vector<size_t> coord(1, 1);
		size_t pole = 0;
		while(pole != 1) {
			pole = 0;
			if(fabs(this->get_F(coord, pre_N)) < _epsilon) {
				is_stop = true;
				break;
			}
			coord[pole]++;
			while(coord[pole] = (1 << pre_N)) {
				coord[pole] = 1;
				pole++;
				if(pole == 1)
					break;
				coord[pole]++;
			}
		}

		this->__U.clean(pre_N);
	}
	if(pre_N == _M)
		this->solve_coarest(pre_N);
	else {
		for(size_t i = 0; i != nu2; i++) {
			this->smooth(pre_N);
		}
	}
	if(_N > _M)
		is_stop = false;
	while(!is_stop) {
		this->putback_e(pre_N);
		pre_N++;
		for(size_t i = 0; i != nu2; i++) {
			this->smooth(pre_N);
		}
		if(pre_N == _N) {
			is_stop = true;
		}
	}
};

#else
//Do nothing.
#endif
